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Abstract. A contact discontinuity tracking method with a specially designed moving grid is developed to 
eliminate the interface smearing completely. In order to precisely locate the contact surface, an updated 
Riemann solver for unsteady one-dimensional inviscid flows is also developed to allow consideration of the 
specific heat ratio change across the shock wave. These two new computational techniques are illustrated 
in a high Mach number shock tube flow field computation. 
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Nomenclature 

C: sound speed 

C p : specific heat at constant pressure 
C v : specific heat at constant volume 
e: internal energy 
e v p. vibrational energy of species i 
h: enthalpy 
p : pressure 

qp. mass fraction of species i 
Rj : gas constant of species i 
t: time 

T: temperature 

u: velocity 

U : general variable 

IV: shock Mach number 

x : coordinate along the shock tube 

Greek symbols 

( Ah{)T ■ standard heat of formation of species i at tem- 
perature T 

net production rate of species i 
7: specific heat ratio, C p /C v 
gp. molecular weight of species i 
9 V i : characteristic vibrational temperature of species i 
g: density 

r v p. characteristic relaxation time of species i 



Superscript 

*: intermediate states 
first order derivative 
second order derivative 

Subscript 

1: left side 
r: right side 

1 Introduction 

Renewed interest in hypersonics has created a demand for 
short duration high-enthalpy ground test facilities which 
are confined to hypervelocity flight regions. This special 
requirement has spawned a number of exotic facilities 
which represent in essence various extensions of the basic 
shock tube flow. The present work investigates shock tube 
and shock tunnel (without piston) flows, which represent 
only a few of many different types of impulse facilities. 
Numerical simulations of such flows are required to bet- 
ter understand the test flow conditions, to determine test 
duration time and supplement the results of experiments. 
It is also useful to improve the shock tunnel operation 
characteristics. 

However, flow conditions in a shock tunnel are very 
severe and put enormous strain on the accuracy and sta- 
bility of the current numerical techniques. The unsteady 
flow inside a high enthalpy shock tube and shock tun- 
nel is significantly influenced by very strong shock waves, 
chemically reacting and inert gas interfaces (contact dis- 
continuities) , high temperature effects, chemical reactions 
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and viscous interactions. These important flow features 
take place in a few millimeters in the axial direction in a 
facility that can be tens of meter long. 

A number of researchers (see, for example, Refs. 1- 
4) have performed numerical simulations of high enthalpy 
pulse facilities to tackle the above mentioned challenges. 
Cambier et al. (1992) describe the computational work 
on the flow simulation of the NASA Ames 16” Shock 
Tunnel Facility and describe the numerical problems en- 
countered during the computation of various flow tran- 
sients and the methods used to resolve them. Particular 
attention is given to problems arising from extremely stiff 
chemistry. Wilson (1992) presents a quasi-one-dimensional 
methodology for numerically simulating the flow inside 
high-entlralpy pulse facilities. The numerical approach 
uses a finite volume TVD scheme for the Euler equations 
coupled with finite- rate chemistry on a moving mesh. A 
Riemann solver is incorporated for tracking contact dis- 
continuities. A 4tlr order non-oscillatory scheme and flow 
loss models are used by Itolr et al (1993) to study the 
tuned operation of a free piston shock tunnel. The aim 
of the present paper was to develop a numerical simula- 
tion technique for shock tube and shock tunnel flows with 
particular emphasis on predicting accurate test times by 
accurately determining the contact surface motion sep- 
arating the chemically active region of the flow from the 
cold flow and taking into account the real gas effects on the 
shock wave strengths with an improved Riemann solver. 

The numerical approach taken in the present work has 
been to solve the unsteady quasi- one-dimensional Euler 
equations coupled to a detailed finite rate chemistry model 
for high temperature air. The numerical scheme employed 
is a modified version of the two-dimensional Godunov type 
upwind non-oscillatory Euler solver proposed by Rodionov 
(Rodionov 1987) for transient flows. Second order spatial 
accuracy of smooth solutions is obtained by linear approx- 
imation of the conservative variables within the control 
volumes. The corresponding gradient of the flow variables 
in the control volume determined by one sided deriva- 
tives is chosen as the minimum of average gradients for 
neighbouring control volumes (Colgan 1972). Second or- 
der accuracy with respect to time is achieved by a two step 
predictor-corrector technique. In the predictor step, inter- 
mediate values at grid points are determined following the 
above mentioned one-sided derivatives in conjunction with 
the principle of minimum derivative. Using the arithmetic 
average of the known variable value and the correspond- 
ing intermediate value to replace the known value of grid 
points an exact Riemann solver is used in the corrector 
step to determine the flow variables at the control volume 





Fig. 2. Four available wave patterns in the Riemann prob- 
lem solution S - shock wave R - rarefaction wave C - contact 
surface 



interfaces. Note that the predictor step of the scheme in- 
volves no Riemann solver and hence contributes to the 
efficiency of the employed computational technique. 

Two numerical issues are the focus of the present work. 
The first issue relates to the solution of the exact Rie- 
mann problem for very high Mach number shock tube or 
shock tunnel flows involving very strong discontinuities 
when the specific heat ratio becomes a function of tem- 
perature. In these cases, it is inconsistent to use an exact 
Riemann solver based on constant specific heat ratio when 
chemical reactions are treated as temperature dependent. 
An improved exact Riemann solver for unsteady one di- 
mensional flows is presented which takes into account the 
change of the specific heat ratio across shock waves. 

The second issue concerns the smearing of cold and 
very hot gas interfaces by numerical diffusion. This smear- 
ing increases as the distance the interface travels becomes 
greater. This is the case when the length of the shock tun- 
nel being simulated becomes large. This large smearing 
results in non-physical mass transfer of chemical species 
and in non-physical chemical reactions across the inter- 
face. Moreover, the performance of the shock tube and 
the available maximum test time in the test section of the 
shock tunnel are very closely related to the interaction 
between the reflected shock wave from the end of of the 
driven section and the contact discontinuity, the interface 
between cold and reacting hot gases. Therefore, it is im- 
perative to precisely locate the contact discontinuity at 
any instant of time in the shock tube. 

Tracking the contact discontinuity eliminates the inter- 
face smearing completely. A contact discontinuity tracking 
method is proposed in the present work by applying the 
above described computational technique to a specially 
designed moving grid, resulting in sharp, step-like varia- 
tions of the flow parameters across it. 
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2 An updated Riemann solver 

Consider the grid shown in Fig. 1 with Uj as the known 
value of the general variable U at grid point Xj, and time 
level t n . In order to find the value of Uj at time level t n+1 , 
Uj 1 , with a finite volume method, a certain discretization 
scheme is needed to determine the values of the variable U 
at the cell boundaries, Uj_\/ 2 and f/ ?+1 / 2 - However, while 
discontinuities exist in the flow held, the discontinuity be- 
tween two adjacent points, say Uj at Xj and Uj + 1 at Xj+i, 
will break into leftward and rightward moving waves sep- 
arated by a contact surface. Each wave can be either a 
shock wave or a rarefaction wave and the physically avail- 
able combinations produce four wave patterns as shown 
in Fig. 2. The problem of determining the wave patterns 
and the how held in each region is called a Riemann prob- 
lem and the algorithm to solve this problem is called a 
Riemann solver. Therefore, while discontinuities exist in 
the how held, in order to determine the variable values 
at the cell boundary, a Riemann problem must be solved. 
Thus an efficient and precise Riemann solver is an impor- 
tant part of a CFD method for supersonic and hypersonic 
hows where the how discontinuities must be considered. 

Riemann problems were hrst introduced into computa- 
tional fluid dynamics by Godunov (Godunov 1976). Since 
then, various approximate solvers and “exact” solvers were 
developed (see for example, Gottlieb, Groth 1988). 

The exact Riemann solver developed by Gottlieb and 
Groth is demonstrated to be an efficient and robust one. 
The central concept of the solver is to choose the common 
how velocity of the intermediate state (tt*) as the iterate 
and the pressure difference (p* — p*) across this contact 
discontinuity is made equal to zero. 

The initial guess (n= 0) for the how velocity u* is given 

by 

(1) 



l + z 



where 



u t =Ul + 



u = u r + 



2 Ci 
7 / - 1 
2 C r 

lr - 1 



z = 7 1 - 1 ( V Pi V 
7r - 1 Cl P r 
if Pi <Pr, 7 = lr 
Pi >Pr, 7 = ll ■ 

By using the standard Newton iteration procedure the 
successive iterates of u* are given by 
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where s p is a prescribed infinitesimal (le-7 in the present 
work). In order to improve the convergence rate, second 
order derivatives are also included in the Newton iterative 
procedure. The successive iterates of u* will then be: 



*n+ 1 
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(Pl(<) ~Pr«)) 2 (Pl" «)-Pr"«)) 
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where the double prime denotes the second derivative with 
respect to u*. 

The detailed expressions for the four wave patterns in 
Fig. 2 are given in [8] (Gottlieb, Groth 1988). These ex- 
pressions are derived based on the assumption that the 
specihc heat ratio 7 is constant across the shock wave. In 
high Mach number shock tube and shock tunnel or any 
other flow held with strong discontinuities, the tempera- 
ture variation in the held is large, and the specihc heat 
ratio 7 can not be considered constant as it is a function 
of temperature. Sometimes, even the chemical reactions 
must be considered. It is inconsistent to use the Riemann 
solver based on constant specihc heat ratio assumption 
while the main program considers chemical reactions and 
treats the thermodynamic properties as a function of tem- 
perature. As an update to the mentioned Riemann solver, 
further iteration is introduced to consider the specihc heat 
ratio change across the shock wave, if the shock wave 
strength in the Riemann problem solution is very strong, 
for example, the pressure ratio across the shock is greater 
than 10 in the present work. 

For non-constant specihc heat ratio, we have 
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otherwise 
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Pi 



where 7 = 



H + it 



where the prime denotes the hrst derivative with respect 
to u*. The values of pf(u*), p*«), pfiK), pf (K) are 
determined based on the wave patterns. The iteration ends 
when the following convergence criterion is satished: 



Fewer mathematical operations will be needed for itera- 
tions to consider the effect of non-constant specihc heat 
ratio, if the common pressure of the intermediate states 
(p*) is chosen as an iterate. The initial guess of pressure 
p* should obviously be 



P* = ( Pl+Pr)/2 



(!') 
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For the guessed p* , based on constant 7 calculate an as- 
sumed T* and then find the guessed 7*. Once the guessed 
p* and 7* are prescribed, the T* can then be recalculated. 
The successive iterates of p* is as follows. 



Pn+1 Pn 



U*(P*n)-<(Pn) 
U T (Pn) ~ U *' (Pn) 



(Ul(Pn) - K(Pn)) 2 (Ui" ( Pn ) ~ U*r" (p*)) 

2{u* l '(p* n )-uV(p* n )Y 



(4') 
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where 



The iteration ends when a certain convergence criterion is 
satisfied: 

I 2« - u* r ) I 



« + <) l<£ “’ (30 
where e u is a prescribed infinitesimal (0.001 in the present 
work). All the corresponding terms involved in the calcu- 
lation for the four wave patterns can be written as follows 
when the specific heat ratio change across the waves is 
considered. 

1. In the case of rightward moving shock for which p* > p r , 
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For the cases of rarefaction waves, the arithmetic average 
specific heat ratio is used and the terms involved in the 
calculation are shown in case 3 and 4. 

3. In the case of right rarefaction wave with p* < p r , 



K = u r + -{C* - C r ), 

7-1 
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4. In the case of left rarefaction wave with p* < pi, 

2 
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2. In the case of leftward moving shock for which p* > p;, 
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I.S - Incident Shock wave 
R - Rarefaction wave 
C - Contact surface 
R.S - Reflected Shock wave 
R.R - Reflected Rarefaction wave 



Fig. 3. Shock tube scheme 




ti n+1 
"b n 

Fig. 4. Grid distribution 



3 Shock tube flow field numerical simulation 

In shock tube flow field analysis, it is important to pre- 
cisely locate the contact surface. A sketch of a shock tube 
is shown in Fig. 3. When the diaphragm breaks, a strong 
incident shock followed by a contact surface is produced 
by the moving high pressure gas in the driver towards the 
lower pressure driven tube. This shock wave travels down 
to the end of the driven tube and then reflects back to 
interact with the contact surface. The available maximum 
test time is closely related to the interaction between the 
reflected shock wave and the contact surface. In order to 
obtain high Mach number and stagnation temperature in 



the test section, hydrogen (or helium) with high pressure 
is used as the driver gas; the driven gas is air (78.084% 
N 2 , 20.946% O 2 , 0.97% A,,) at standard atmosphere pres- 
sure. In such a case, the incident shock may be very strong, 
and the temperature behind the shock is high enough that 
the chemical reactions and non-constant thermal property 
effect must be considered. A Riemann solver based on 
constant specific heat ratio assumption is not acceptable. 
In order to illustrate the appropriateness of the updated 
Riemann solver, a high Mach number shock tube flow is 
numerically solved. The shock tube scheme and detailed 
input data are shown in Fig. 3. 
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Density Contour (log D) Pressure Contour (log P) 




Mach Number Contour Temperature Contour (log T) 




3.1 Governing differential equations 



Since the incident shock wave is quite strong, chemical 
reactions must be considered for the driven gas. In this 
work, five species in air (O 2 , N 2 , NO, N, O) are involved 
and seventeen chemical reactions are considered. The re- 
lated thermodynamic data and chemical reaction infor- 
mation are taken from Ref. 9 and Ref. 10. The governing 
differential equations for a non-equilibrium unsteady one 
dimensional flow in a shock tube problem can be written 
as follows: 



dg dgu 
dt + dx 



( 5 ) 



d(gu) d(p+gU 2 ) 



= 0, 



d u 2 

-Q t [Q{h + Zne V i + y ) - P\ 
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0 d 

g- t mi) + 

i^(mevi) + -T^iguqievi) = Qe V i)$i 



+gq, 



[0 vi Ri]/[ex p(%4) - 1] - e v 



where 
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Hi 2 3 4 
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~[Ah{] 0 ) 



dt 



dx 
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Table 1. Chemical reactions involved 





A(i,m)X(i) 


+ A(i,m)X(i) 


== B(i,m)X(i) 


+B(i,m)x(i) 


+B(i,m)X(i) 


1 


o 2 


+ N 


== 20 


+ N 




2 


o 2 


+ NO 


== 20 


+ NO 




3 


o 2 


+ o 2 


== 20 


+ o 2 




4 


o 2 


+ n 2 


== 20 


+ n 2 




5 


o 2 


+ o 


== 20 


+ o 




6 


n 2 


+ N 


== 2N 


+ N 




7 


n 2 


+ NO 


== 2N 


+ NO 




8 


n 2 


+ o 2 


== 2N 


+ o 2 




9 


n 2 


+ n 2 


== 2N 


+ N 2 




10 


n 2 


+ o 


== 2NO 


+ o 




11 


NO 


+ n 2 


== N 


+ o 


+ n 2 


12 


NO 


+ N 


== N 


+ o 


+ N 


13 


NO 


+ o 2 


== N 


+ o 


+ o 2 


14 


NO 


+ o 


== N 


+ o 


+ 0 


15 


NO 


+ NO 


:== N 


+ o 


+ NO 


16 


0 


+ NO 


== N 


+ o 2 




17 


0 


+ n 2 


== N 


+ NO 





C p i — 8314.0(aii + a^T + a^T 2 + a^T 3 + a^T 4 ), 

-A(i, m)){AFR{m)T BFR{rn) 
cxp( DFR ^> ) JJ t ( 1000 ^ \ARm) 

- ABR(m)T BBR WeM - DB ^ ) 

jj / 1000/b^j . l 

£> 

^^exp(A 4 [T-i - 0.015(^)i] - 18.42) 

T vi p y qp_ 1 

1.033e5^ fc Atfe 

Constants in the vibrational energy equations are as fol- 
lows: 



Species 


n 2 


o 2 


NO 


A-i 


220 


129 


168 


@vi 


3390 


2270 


2740 



Constants related to chemical reactions see tables and 
other data can be found in [9] (Sagnier, Marraffa 1991) 
and [10] (Gupta, Yos, Thompson, Lee 1990). 

3.2 Contact surface tracking 

In a shock tube or a shock tunnel, the driver and the driven 
gases are usually different. When the diaphragm breaks, a 
contact surface separates the two different gases. Tracking 
the contact surface to divide the computational domain 
into two parts with two different gas species can greatly 
simplify the computation work. It is especially important 
when the incident shock is strong and the temperature is 
high enough behind the incident shock wave that the real 
gas effect must be considered. When the incident shock 
reaches the tube end and a reflected shock wave is devel- 
oped, the whole driven gas is under high temperature, the 



Table 2. Chemical reaction rate coefficients 





AFR 


BFR 


DFR 


ABR 


BBR 


DBR 


1 


3.6el8 


-1 


0.595e5 


3.0el5 


-0.5 


0.0 


2 


3.6el8 


-1 


0.595e5 


3.0el5 


-0.5 


0.0 


3 


3.249el9 


-1 


0.595e5 


2.7el6 


-0.5 


0.0 


4 


7.2el8 


-1 


0.595e5 


6.0el5 


-0.5 


0.0 


5 


9.0el9 


-1 


0.595e5 


7.5el6 


-0.5 


0.0 


6 


4.08e22 


-1.5 


1.132e5 


2.27e21 


-1.5 


0.0 


7 


1.9el7 


-0.5 


1.132e5 


l.lel6 


-0.5 


0.0 


8 


1.9el7 


-0.5 


1.132e5 


l.lel6 


-0.5 


0.0 


9 


4.7el7 


-0.5 


1.132e5 


2.72el6 


-0.5 


0.0 


10 


1.92el7 


-0.5 


1.132e5 


l.lel6 


-0.5 


0.0 


11 


3.97e20 


-1.5 


0.755e5 


1.0e20 


-1.5 


0.0 


12 


7.8e20 


-1.5 


0.755e5 


2.0e20 


-0.5 


0.0 


13 


3.97e20 


-1.5 


0.755e5 


1.0e20 


-1.5 


0.0 


14 


7.8e20 


-1.5 


0.755e5 


2.0e20 


-1.5 


0.0 


15 


7.8e20 


-1.5 


0.755e5 


2.0e20 


-1.5 


0.0 


16 


3.18e9 


1.0 


1.968e4 


1.3el0 


1.0 


3.58e3 


17 


7.0el3 


0.0 


3.8e4 


1.56el3 


0.0 


0.0 



real gas effects and the chemical reactions must be con- 
sidered for the whole region. In the meantime, the driver 
gas is under low temperature, even though the starting 
temperature may be higher than the one in driven gas be- 
cause of the rarefaction wave caused by the breakdown of 
the diaphragm. The chemical reactions there can be ne- 
glected. Tracking the contact surface to separate these two 
regions can limit the consideration of real gas effects to the 
driven gas only. In addition, the performance of the shock 
tube and the available maximum test time, are closely re- 
lated to the interaction between the reflected shock wave 
and the contact surface. Therefore, it is also important to 
precisely locate the contact surface. 

Obviously, at time level t°= 0, the starting point of 
the contact surface is at the location of the diaphragm, 
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Fig. 6a, b. Field parameter variations at different time levels a current method, b conventional method 
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Fig. 6a, b. (continued) 
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x = Xj S . |_i /2 = Xj S + Ax/2 as shown in Fig. 4. On the left 
side is the uniform driver gas and on the right side is the 
uniform driven gas; solving the Riemann problem gives 
the location of the contact surface, the angle, a, in the 
Fig. 4. Once the time interval, r, is set, the starting point 
for the new time level t 1 = t° + r is determined. At this 
new time level, the left-hand side and the right-hand side 
of the contact surface are not uniform any longer. The ini- 
tial values of the Riemann problem with no consideration 
of gradients on both sides of the discontinuity are deter- 
mined by the discretization scheme. Solving the Riemann 
problem will give the location of the contact surface (angle 
a) at the new time level. In this way, the contact surface 
can be tracked for all time levels. The precision of the 
contact surface depends on the precision of the Riemann 
solver and the discretization scheme. 



3.3 Grid distribution and discretization 



The whole solution domain is divided into cells as shown 
in Fig. 4. 

For each cell, the differential equations can be con- 
verted to algebraic equations using Green’s theorem. Ar- 
ranging the grid line distribution by tracking the contact 
surface can greatly improve the computational precision. 
Obviously, at time level t n = 0, the starting point of 
the contact surface is at the location of the diaphragm, 
x = Xj S+ 1/2 = Xj S + Ax/2 as shown in Fig. 4. The grid 
intervals in the driver section and the driven section are 
equal, Axl = Axr = Ax. According to the Godunov- 
Colgan method (Colgan 1972), an approximation using 
one-sided derivatives in conjunction with the principle of 
the minimum derivative provides a correct distribution of 
variables in the cells when discontinuities are present, and 
is of 0 (h 2 ) accurate for smooth solutions. Based on the 
principle of the minimum derivative, the independent vari- 
able values U at the cell boundary js + 1/2 from the left 
cell ,js, should be 



U js- t-i/2 — U j s A j s * Axl/ 2; 



Kjs 



min( 



U'js + l U'js 




v* js U'js— 1 


Xjs-\- 1 Xjs 


7 


Xjs Xjs— 1 



)• 



The independent variable values U at the cell boundary 
js + 1/2 from the right cell, js + 1, should be 



Ujs- j-1/2 — Ujs+ 1 Ay s _|_i * Ax Ft/ 2, 



Kjs+i = min( 



^js-\- 2 'Ujs+l 




U'js + l V* js 


Xjs-\- 2 ^js + 1 


7 


Xjs-\- 1 Xjs 



)• 



Thus, a discontinuity exists at the cell boundary. A Rie- 
mann problem based on these initial values with no con- 
sideration of gradients on both sides of the discontinuity 
is solved to determine the contact surface location angle, 
a , and the cell boundary value. The maximum velocity 
of propagation of disturbances generated in the Riemann 
problem, C w , is also recorded for the purpose of deter- 
mining the time interval, r. This discretization scheme is 



proved to be first order accurate in marching direction 
(time in this work) and second order accurate in the coor- 
dinate x. It is stable if the following criterion is satisfied. 



t/Ax <= 0.5 /C w , 

where C w is the maximum velocity of propagation of dis- 
turbances generated in the Riemann problem and r is the 
permissible time interval. Once the contact surface loca- 
tion angle, a, and the time interval, r, are determined, the 
grid distribution for the new time level t n+1 = t n + t can 
be found geometrically as shown in Fig. 4. 

In order to improve the accuracy in the marching direc- 
tion, a modification developed by Rodionov is introduced 
(Rodionov 1987). According to Rodionov, intermediate 
values at grid points are determined following the above 
mentioned one-sided derivatives in conjunction with the 
principle of the minimum derivative before solving the Rie- 
mann problem. Then the arithmetic average of the known 
variable value and the corresponding intermediate value 
is used to replace the known value at grid points and to 
discretize the equations. This new scheme can provide an 
accuracy close to second order in the marching direction. 
When the grid distribution is determined, one follows the 
same procedure to solve the Riemann problems for all the 
cell boundaries to determine the cell boundary values and 
the maximum velocity of propagation of disturbances gen- 
erated in the Riemann problems, C w . The final discretized 
equations can be written as follows: 

(6<) 

{eu)j =-( ( Qu)j + J- ( [p + Qu(u - W*)\ 1 

c \ ll x \ J 2 /g/\ 

-[p+ gu(u-W*)\ j+ ^, 
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C = ( W j+l/2 - W*_ 1/2 )T/h x , 

5 = 0.5, 

51 = 0.5, 

5 2 = 0.4. 



3.4 Computational results 

The high Mach number shock tube problem outlined in 
Fig. 3 is solved with the contact discontinuity tracking 
method and the updated Riemann solver. The results are 
summarized in Fig. 5 and Fig. 6. Figure 5 shows the con- 
tour picture of pressure, temperature, density and Mach 
number. In order to assess the quantitative variations in 
flow variables along the shock tube, the parameter varia- 
tions at three typical time levels are shown in Fig. 6A. The 
time level t= 0.5 msec, corresponds to the stage when the 
incident shock and the rarefaction waves have not reached 
the tube end; Time level t= 1 msec, corresponds to the 
stage after the rarefaction waves have reached the tube 



end, but the incident shock has not. Time level t= 1.6 
msec, displays the case after the incident shock reaches 
the wall and reflects from there. The same shock tube 
problem is also solved by a conventional method with no 
contact surface tracking and no chemical reaction consid- 
ered. The specific heat ratio is considered to be constant 
both in the main flow field computation and in the Rie- 
mann solver. The results are shown in Fig. 6B. Comparing 
these results, it can be seen that both methods can pre- 
dict the incident shock by the sharp jump of pressure and 
density profiles. However, only the present method can lo- 
cate the contact surface by the sharp jump of the density 
profile. There is a large smeared region around the contact 
surface which makes it difficult to locate the contact sur- 
face when the conventional method is used. Since the real 
gas effects are considered in the current method, there is 
no sharp jump for the reflected shock wave. The reflected 
shock wave travels in a non-uniform reacting fluid inter- 
medium, no sharp jump across the reflected shock will 
appear. The sharp jump shown in Fig. 6B is caused by the 
perfect gas assumption, it does not mean the conventional 
method can precisely locate the reflected shock wave. For 
the field behavior, the shock wave strength and location, 
there are significant differences between the predicted re- 
sults by these two methods. This is caused by real gas 
effects considered in the present method. The computa- 
tional technique such as the contact discontinuity track- 
ing method and the Riemann solver can change the solu- 
tion precision and the computation cost; however, only the 
equations, including the chemistry models, viscous mod- 
els will influence the field behavior. Therefore, the detailed 
discussion of the field behavior is beyond the scope of this 
article. 



4 Conclusions 

The Riemann solver developed by Gottlieb and Groth in 
1987 is demonstrated to be an efficient and robust one. 
Since it is based on constant specific heat ratio assump- 
tion across the shock wave, it is not acceptable when a 
strong discontinuity is present in the flow field. Espe- 
cially, when chemical reactions are considered and ther- 
modynamic properties are treated as functions of temper- 
ature in the main code, it is inconsistent to take the con- 
stant specific heat ratio assumption in the Riemann solver 
subroutine. An update to consider the specific heat ratio 
change across shocks in the solver must be introduced. 
However, even if there are strong multi-discontinuities in 
the whole flow field, most calls for the Riemann solver in 
the code only handle a weak shock in the Riemann prob- 
lem, therefore, the update is introduced in the solver as 
an option only when it is necessary. This treatment can 
keep the merit of the original solver and simplify the mod- 
ification. The solver starts with the constant specific heat 
ratio assumption and then processes further iterations to 
consider the effect of variable specific heat ratio across the 
shock wave when it is necessary. In this part of further it- 
erations, the conventional way of choosing the common 
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pressure of the intermediate states (p*) as the iterate is 
more convenient. 

In shock tube or shock tunnel, the smearing of cold and 
hot gas interfaces by numerical diffusion increases as the 
distance that the interface travels becomes greater. Track- 
ing the contact discontinuity can eliminate the interface 
smearing completely so that the simulation precision is 
greatly improved. Besides, in high Mach number shock 
tube or shock tunnel, the driver gas and the driven gas 
are usually different. When the diaphragm breaks, a con- 
tact surface separates the two different gases. Tracking the 
contact surface to divide the computational domain into 
two parts with two different gas species can greatly sim- 
plify the computational work. It is especially important 
when the incident shock is strong so the real gas effects 
and the chemical reactions must be considered. Tracking 
the contact surface to separate these two regions can limit 
the consideration of real gas effect in the driven gas only. 
Obviously, similar tracking methods based on the perfect 
gas assumption [11] [12] do not provide this advantage and 
their applications are therefore limited. 



References 

Cambier JL, Tokarcik S, Prabhu DK (1992) Numerical simula- 
tions of unsteady flow in a hypersonic shock tunnel facility. 
AIAA Paper 92:4029 

Wilson G (1992) Time-dependent quasi-one-dimensional sim- 
ulations of high enthalpy pulse facilities. AIAA Paper 
92:5096 



Itoh K, Tani K, Tanno H, Takahashi M, Miyajima H, Asano 
T, Sasoh A, Takayama K (1993) A numerical and experi- 
mental study of the free piston shock tunnel. Proc 19th Int 
Symp Shock Waves, pp 257-263 

Jacobs PA (1991) Simulation of transient flow in a shock tunnel 
and a high Mach number nozzle. Proc 4th Int Symp Fluid 
Dynamics 

Rodionov AV (1987) Second order accurate monotonic scheme 
for nonequilibrium flow computation. J of Comp Math and 
Math Phys 27 

Colgan VP (1972) Application of the principle of minimum val- 
ues of the derivative to the construction of finite difference 
schemes for the discontinuous solution of gasdynamics. Sci- 
entific Notes, TsAGI, No 6 

Godunov SK (1976) Numerical solution of multidimensional 
problems in gasdynamics. Nauka, Moscow 

Gottlieb J J, Groth PT (1988) Assessment of Riemann solvers 
for unsteady one-dimensional inviscid flows of perfect gases. 
J of Comp Phys Oct, pp 437-458 

Sagnier P, Marraffa L (1991) Parametric study of thermal 
and chemical nonequilibrium nozzle flow. AIAA J March, 
pp 334-343 

Gupta RN, Yos JM, Thompson RA, Lee K (1990) A review of 
reaction rates and thermodynamic and transport proper- 
ties for an 11-species air model for chemical and thermal 
nonequilibrium calculations to 30000K. NASA Reference 
Publication 1232 

Jacobs PA (1994) Quasi-one-dimensional modeling of a free- 
piston shock tunnel. AIAA J, pp 137-145 

Wilson D, Tan Z, Varghese P (1996) Numerical simulation of 
the blast-wave accelerator. AIAA J, pp 1341 1347 




